Path Integral Monte Carlo Calculation of the Deuterium Hugoniot 
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Restricted path integral Monte Carlo simulations have been used to calculate the equilibrium 
properties of deuterium for two densities: 0.674 and 0.838 gcm""^ (ts = 2.00 and 1.86) in the 
temperature range of 10 000 K < T < 1 000 000 K. Using the calculated internal energies and 
pressures we estimate the shock hugoniot and compare with recent Laser shock wave experiments. 
We study finite size effects and the dependence on the time step of the path integral. Further, we 
compare the results obtained with a free particle nodal restriction with those from a self-consistent 
variational principle, which includes interactions and bound states. 



PACS Numbers: 71.10.-w 05.30.-d 02.70.Lq 



I. INTRODUCTION 



Recent laser shock wave experiments on pre- 
compressed liquid deuterium ||l|,|^ have produced an un- 
expected equation of state for pressures up to 3.4 Mbar. 
It was found that deuterium has a significantly higher 
compressibility than predicted by the semi-empirical 
equation of state based on plasma many-body theory 
and lower pressure shock data (see SESAME model [^). 
These experiments have triggered theoretical efforts to 
understand the state of compressed hydrogen in this 
range of density and temperature, made difficult be- 
cause the experiments are in regime where strong cor- 
relations and a significant degree of electron degeneracy 
are present. At this high density, it is problematic even 
to define the basic units such as molecules, atoms, free 
deuterons and electrons. Conductivity measurements ||] 
as well as theoretical estimates suggest that in the 
experiment, a state of significant but not complete met- 
alization was reached. 

A variety of simulation techniques and analytical mod- 
els have been advanced to describe hydrogen in this par- 
ticular regime. There are ab initio methods such as re- 
stricted path integral Monte Carlo simulations (PIMC) 
and density functional theory molecular dynamics 
(DFT-MD) |,|. Further there are models that min- 
imize an approximate free energy function constructed 
from known theoretical limits with respect to the chemi- 
cal composition, which work very well in certain regimes. 
The most widely used include | p^ , pl] ,p[ . 

We present new results from PIMC simulations. What 
emerges is a relative consensus of theoretical calculations. 
First, we performed a finite size and time step study us- 
ing a parallelized PIMC code that allowed simulation of 
systems with Np — 64 pairs of electrons and deuterons 
and more importantly to decrease the time step from 



II. RESTRICTED PATH INTEGRALS 

The density matrix of a quantum system at temper- 
ature ksT =1/(3 can be written as a integral over all 
paths Rt, 

Rt stands for the entire paths of N particles in 3 dimen- 
sional space Rt = {ru, ■ ■ ■ ,i^Nt) beginning at Rq and 
connecting to VTip. V labels the permutation of the par- 
ticles. The upper sign corresponds to a system of bosons 
and the lower one to fermions. For non-relativistic parti- 
cles interacting with a potential ^(R), the action of the 
path S\R.t\ is given by. 
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8 ■ 10 K. More importantly, we 



studied the effect of the nodal restriction on the hugoniot. 



One can estimate quantum mechanical expectation val- 
ues using Monte Carlo simulations with a finite num- 
ber of imaginary time slices M corresponding to a time 
step T = 13/M. 

For fermionic systems the integration is complicated 
due to the cancellation of positive and negative contribu- 
tions to the integral, {the fermion sign problem). It can 
be shown that the efficiency of the straightforward imple- 
mentation scales like e~'^^^^ , where / is the free energy 
difference per particle of a corresponding fermionic and 
bosonic system [|l3|. In |14 1^, it has been shown that 
one can evaluate the path integral by restricting the path 
to only specific positive contributions. One introduces a 
reference point R* on the path that specifies the nodes 
of the density matrix, p(R, R*,t) = 0. A node-avoiding 
path for < i < /3 neither touches nor crosses a node: 
p(R(t),R*,i) ^ 0. By restricting the integral to node- 
avoiding paths. 
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Jd-Ro pf(Ro,R*;0) ^ dRt e-^[^*l, (3) 



R.0 — GT(R,* 



(T(R*) denotes the restriction) the contributions are 
positive and therefore PIMC represents, in principle, a 
solution to the sign problem. The method is exact if the 
exact fermionic density matrix is used for the restriction. 
However, the exact density matrix is only known in a 
few cases. In practice, applications have approximated 
the fermionic density matrix, by a determinant of single 
particle density matrices. 



p(R,R';/3) = 



pi(ri,ri;/3) ... pi{rN,r[;/3) 
pi(ri,r^;/3) ... pi(rAr, r'^; /3) 



(4) 



This approach has been extensively applied using the free 
particle nodes [p^ , 

Pi(r,r',/3) = (47rA/3)-3/2exp{-(r-r')V4A/3} (5) 

with A = ?i^/2m, including applications to dense hydro- 
gen . It can be shown that for temperatures larger 
than the Fermi energy the interacting nodal surface ap- 
proaches the free particle (FP) nodal surface. In addi- 
tion, in the limit of low density, exchange effects are negli- 
gible, the nodal constraint has a small effect on the path 
and therefore its precise shape is not important. The 
FP nodes also become exact in the limit of high density 
when kinetic effects dominate over the interaction poten- 
tial. However, for the densities and temperatures under 
consideration, interactions could have a significant effect 
on the fermionic density matrix. 

To gain some quantitative estimate of the possible ef- 
fect of the nodal restriction on the thermodynamic prop- 
erties, it is necessary to try an alternative. In addition to 
FP nodes, we used a restriction taken from a variational 
density matrix (VDM) that already includes interactions 
and atomic and molecular bound states. 

The VDM is a variational solution of the Bloch equa- 
tion. Assume a trial density matrix with parameters qi 
that depend on imaginary time /3 and R', 



p(R,R';/3) =p(R,gi,. 



(6) 



By minimizing the integral: 
/ap(R,R';/3) 



dR 



V d(3 



+ np{R,R';p) =0 , (7) 



one determines equations for the dynamics of the param- 
eters in imaginary time: 

IdH ^ ■ r 

■7;^+ M q^Q where H= pHp dR . (8) 

z oq J 

The normalization matrix is: 



dqidq'. 



JdRp{R,q;P) p{R,q';P) 



(9) 



We assume the density matrix is a Slater determinant of 
single particle Gaussian functions 

pi(r,r',/3) = (ttw)-^/^ exp {-{r - mf /w + d} (10) 

where the variational parameters are the mean m, 
squared width w and amplitude d. The differential equa- 
tion for this ansatz are given in |p^ . The initial condi- 
tions at /? — > are w = 2/3, m = r' and c? = in order 
to regain the correct FP limit. It follows from Eq. 
that at low temperature, the VDM goes to the lowest 
energy wave function within the variational basis. For 
an isolated atom or molecule this will be a bound state, 
in contrast to the delocalized state of the FP density 
matrix. A further discussion of the VDM properties is 
given in . Note that this discussion concerns only the 
nodal restriction. In performing the PIMC simulation, 
the complete potential between the interacting charges is 
taken into account as discussed in detail in . 
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FIG. 1. Difference in the internal energy from PIMC sim- 
ulations with VDM and FP nodes vs. temperature using 
Np = 32 and r^^ = 2 ■ lO'^K. 



Simulations with VDM nodes lead to lower internal 
energies than those with FP nodes as shown in Fig. 
Since the free energy F is the integral of the internal 
energy over temperature, one can conclude that VDM 
nodes yield to a smaller F and hence, are the more ap- 
propriate nodal surface. 

For the two densities considered here, the state of deu- 
terium goes from a plasma of strongly interacting but 
un-bound dcuterons and electrons at high T to a regime 
at low T, which is characterized by a significant elec- 
tronic degeneracy and bound states. Also at decreasing 
T, one finds an increasing number of electrons involved in 
long permutation cycles. Additionally for T < 15 625 K, 
molecular formation is observed. Comparing FP and 
VDM nodes, one finds that VDM predicts a higher molec- 
ular fraction and fewer permutations hinting to more lo- 
calized electrons. 
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III. SHOCK HUGONIOT 

The recent experiments measured the shock velocity, 
propagating through a sample of pre-compressed liquid 
deuterium characterized by an initial state, (Eq, Vq, po) 
with T = 19.6 K and po ~ 0.171 g/cm"^. Assuming an 
ideal planar shock front, the variables of the shocked ma- 
terial (E, V, p) satisfy the hugoniot relation p6[ |, 

H = E-Eo + \{V -Vo){p + Po)=(^ . (11) 

We set i?o to its exact value of — 15.886eV per atom |l^ 
and Po — 0. Using the simulation results for p and -E, 
we calculate H{T, p) and then interpolate H linearly at 
constant T between the two densities corresponding to 
rs = 1.86 and 2 to obtain a point on the hugoniot in the 
{p,p) plane. (Results at = 1.93 confirm the function 
is linear within the statistical errors). The PIMC data 
for p, E, and the hugoniot are given in Tab. ||. 
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FIG. 2. Comparison hugoniot function calculated with 
PIMC simulations of different accuracy: FP nodes with 
Np^32 (A for = lO^K reported in [5], [> for 

= 2- lO^K, V for = 8 • lO'^K and = 2- lO^K) and 
Np=6A (□ for r"^ = 2- lO'^K) as well as with VDM nodes and 
Np=32 (o for = lO'^K and • for r"^ = 2 ■ lO^K). Begin- 
ning at high pressures, the points on each hugoniot correspond 
to the following temperatures 125 000, 62 500, 31 250, 15 625, 
and 10 000 K. The dashed line corresponds to a calculation 
using the VDM alone. 

In Fig. we compare the effects of different approx- 
imations made in the PIMC simulations such as time 
step T, number of pairs Np and the type of nodal re- 
striction. For pressures above 3 Mbar, all these approx- 
imations have a very small effect. The reason is that 
PIMC simulation become increasingly accurate as tem- 
perature increases. The first noticeable difference occurs 
at p « 2.7Mbar, which corresponds to T = 62 500 K. At 
lower pressures, the differences become more and more 



pronounced. We have performed simulations with free 
particle nodes and Np = 32 for three different values of 
T. Using a smaller time step makes the simulations com- 
putationally more demanding and it shifts the hugoniot 
curves to lower densities. These differences come mainly 
from enforcing the nodal surfaces more accurately, which 
seems to be more relevant than the simultaneous im- 
provements in the accuracy of the action S, that is the 
time step is constrained more by the Fermi statistics than 
it is by the potential energy. We improved the efficiency 
of the algorithm by using a smaller time step tf for eval- 
uating the Fermi action than the time step tb used for 
the potential action. Unless specified otherwise, we used 
Tp = Tb = T. At even lower pressures not shown in Fig. 
^, all of the hugoniot curves with FP nodes turn around 
and go to low densities as expected. 

As a next step, we replaced the FP nodes by VDM 
nodes. Those results show that the form of the nodes 
has a significant effect for p below 2 Mbar. Using a 
smaller r also shifts the curve to slightly lower densities. 
In the region where atoms and molecules are forming, it 
is plausible that VDM nodes are more accurate than free 
nodes because they can describe those states |l^]. We 
also show a hugoniot derived on the basis of the VDM 
alone (dashed line). These results are quite reasonable 
considering the approximations (Hartree-Fock) made in 
that calculation. Therefore, we consider the PIMC simu- 
lation with the smallest time step using VDM nodes (•) 
to be our most reliable hugoniot. Going to bigger system 
sizes Np = 64 and using FP nodes also shows a shift 
towards lower densities. 
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FIG. 3. Comparison of experimental and several theoreti- 
cal Hugoniot functions. The PIMC curve was calculated with 
VDM nodes, r"^ = 2 ■ lO'' K, and 32 pairs of electrons and 
deuterons. 

Fig. ^ compares the Hugoniot from Laser shock wave 
experiments with PIMC simulation (VDM nodes. 
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= 2 • 10^ K) and several theoretical approaches: 
SESAME model by Kerley § (thin solid line), linear 
model by Ross (dashed line) 0, DFT-MD by 
Lenosky et al. ^ (dash-dotted line), Pade approxima- 
tion in the chemical picture (PACH) by Ebeling et al. 
| pl| (dotted line), and the work by Saumon et al. 
(thin dash-dotted line). 

The differences of the various PIMC curves in Fig. |^ 
are small compared to the deviation from the experimen- 
tal results 0J^. There, an increased compressibility with 
a maximum value of 6 ± 1 was found while PIMC pre- 
dicts 4.3 ±0.1, only slightly higher than that given by the 
SESAME model. Only for p > 2.5Mbar, does our hugo- 
niot lie within experimental errorbars. In this regime, 
the deviations in the PIMC and PACH hugoniot are rel- 
atively small, less than 0.05 gcm~^ in density. In the 
high pressure limit, the hugoniot goes to the FP limit 
of 4-fold compression. This trend is also present in the 
experimental findings. For pressures below 1 Mbar, the 
PIMC hugoniot goes back lower densities and shows the 
expected tendency towards the experimental values from 
earlier gas gun work ||l^,|l^ and lowest data points from 
For these low pressures, differences between PIMC 
and DFT-MD are also relatively small. 



IV. CONCLUSIONS 

We reported results from PIMC simulations and per- 
formed a finite size and time step study. Special emphasis 
was put on improving the fermion nodes where we pre- 
sented the first PIMC results with variational instead of 
FP nodes. We find a slightly increased compressibility 
of 4.3 ± 0.1 compared to the SESAME model but we 
cannot reproduce the experimental findings of values of 
about 6 ± 1 . Further theoretical and experimental work 
will be needed to resolve this discrepancy. 
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TABLE I. Pressure p and internal energy per atom E from PIMC simulations with 32 pairs of electrons and deuterons. For 
T > 250 000 K, we list results from simulations with FP nodes and r^^ = 8 ■ lO** K and r^^ = 2 • 10^ K, otherwise with VDM 



nodes and 


= 2 - lO^'K. 












T{K) 


p{Mbar),rs = 2 


£(eV),r, = 2 


p(Mbar),rs = 1.86 


E{eY),rs = 1.86 


p^"«(gcm-=*) 


p«"9(Mbar) 


1000000 


53.79 ± 0.05 


245.7 ± 0.3 


66.85 ± 0.08 


245.3 ± 0.4 


0.7019 ± 0.0001 


56.08 ± 0.05 


500000 


25.98 ± 0.04 


113.2 ± 0.2 


32.13 ± 0.05 


111.9 ± 0.2 


0.7130 ± 0.0001 


27.48 ± 0.04 


250000 


12.12 ± 0.03 


45.7 ± 0.2 


14.91 ± 0.03 


44.3 ± 0.2 


0.7242 ± 0.0001 


12.99 ± 0.02 


125000 


5.29 ± 0.04 


11.5 ± 0.2 


6.66 ± 0.02 


11.0 ± 0.1 


0.7300 ± 0.0003 


5.76 ± 0.02 


62500 


2.28 ± 0.04 


-3.8 ± 0.2 


2.99 ± 0.04 


-3.8 ± 0.2 


0.733 ± 0.001 


2.54 ± 0.03 


31250 


1.11 ± 0.06 


-9.9 ± 0.3 


1.58 ± 0.07 


-9.7 ± 0.3 


0.733 ± 0.003 


1.28 ± 0.05 


15625 


0.54 ± 0.05 


-12.9 ± 0.3 


1.01 ± 0.05 


-12.0 ± 0.2 


0.721 ± 0.004 


0.68 ± 0.04 


10000 


0.47 ± 0.05 


-13.6 ± 0.3 


0.80 ± 0.08 


-13.2 ± 0.4 


0.690 ± 0.007 


0.51 ± 0.05 
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